
/*!
 * @file mathutil.h
 * @ingroup MathUtil
 *
 * Routines for generating random numbers.
 *
 *
 */

#ifndef _QT_MATHUTIL_H_
#define _QT_MATHUTIL_H_

#include<cstdlib>
#include<cmath>

// template for comments
/*!
 * Description
 *
 * @param
 * @returns
 *
 * History:
 *
 */


namespace qt {

namespace mathutil {


/*!
 * The erfcc function.
 * Adapted adapted from numerical recipies in C.
 * Computes \f[ 2/\sqrt(pi)\int_x^\infty e^(-t^2) dt \f]
 * with fractional error less than 1.2e-7.
 *
 * @param x a positive real number.
 * @return see formula in description.
*/

double erfcc(double x) {

    const double ERFCC_1 = -1.26551223;
    const double ERFCC_2 =  1.00002368;
    const double ERFCC_3 =   0.37409196;
    const double ERFCC_4 =   0.09678418;
    const double ERFCC_5 =  -0.18628806;
    const double ERFCC_6 =   0.27886807;
    const double ERFCC_7 =  -1.13520398;
    const double ERFCC_8 =   1.48851587;
    const double ERFCC_9 =  -0.82215223;
    const double ERFCC_10 =  0.17087277;


    double t,z,ans;
    z=std::abs(x);
    t=1.0/(1.0+0.5*z);
    ans=t*exp(-z*z+ERFCC_1+t*
              (ERFCC_2+t*
               (ERFCC_3+t*
                (ERFCC_4+t*
                 (ERFCC_5+t*
                  (ERFCC_6+t*
                   (ERFCC_7+t*
                    (ERFCC_8+t*
                     (ERFCC_9+t*
                      ERFCC_10)))))))));
    return  x >= 0.0 ? ans : 2.0-ans;
}


/*!
 * The inverse of the erf  function.
 *
 * @param x a real number between -1.0 and 1.0.
 * @return returns the inverse of the erf function, a real number.
*/
double erfinv(double P)
{
    double      Y,A,B,X,Z,W,WI,SN,SD,F,Z2,SIGMA;
    const double      A1=-.5751703,A2=-1.896513,A3=-.5496261E-1;
    const double      B0=-.1137730,B1=-3.293474,B2=-2.374996,B3=-1.187515;
    const double      C0=-.1146666,C1=-.1314774,C2=-.2368201,C3=.5073975e-1;
    const double      D0=-44.27977,D1=21.98546,D2=-7.586103;
    const double      E0=-.5668422E-1,E1=.3937021,E2=-.3166501,E3=.6208963E-1;
    const double      F0=-6.266786,F1=4.666263,F2=-2.962883;
    const double      G0=.1851159E-3,G1=-.2028152E-2,G2=-.1498384,G3=.1078639E-1;
    const double      H0=.9952975E-1,H1=.5211733,H2=-.6888301E-1;
    //  double  RINFM=1.7014E+38;


    X=P;
    SIGMA=(X<0 ? -1.0 : 1.0);

    Z=std::abs(X);
    if (Z>.85) {
        A=1-Z;
        B=Z;
        W=std::sqrt(-1*std::log(A+A*B));
        if (W>=2.5) {
            if (W>=4.) {
                WI=1./W;
                SN=((G3*WI+G2)*WI+G1)*WI;
                SD=((WI+H2)*WI+H1)*WI+H0;
                F=W+W*(G0+SN/SD);
            } else {
                SN=((E3*W+E2)*W+E1)*W;
                SD=((W+F2)*W+F1)*W+F0;
                F=W+W*(E0+SN/SD);
            }
        } else {
            SN=((C3*W+C2)*W+C1)*W;
            SD=((W+D2)*W+D1)*W+D0;
            F=W+W*(C0+SN/SD);
        }
    } else {
        Z2=Z*Z;
        F=Z+Z*(B0+A1*Z2/(B1+Z2+A2/(B2+Z2+A3/(B3+Z2))));
    }
    Y=SIGMA*F;
    return Y;
}


}; // namespace qt::math

}; // namespace qt


#endif
